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ABSTRACT 


A multigrid algorithm is combined with an upwind scheme for solving the two- 
dimensional Reynolds-averaged Navier-Stokes equations on triangular meshes resulting 
in an efficient, accurate code for solving complex flows around multiple bodies. The 
relaxation scheme uses a backward-Euler time difference and relaxes the resulting linear 
system using a red-black procedure. Roe’s flux-splitting scheme is used to discretize 
convective and pressure terms, while a central difference is used for the diffusive terms. 
The multigrid scheme is demonstrated for several flows around single and multielement 
airfoils, including inviscid, laminar and turbulent flows. The results show an appreciable 
speedup of the scheme for inviscid and laminar flows, and dramatic increases in efficiency 
for turbulent cases, especially those on increasingly refined grids. 
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I. INTRODUCTION 


The recent increase in international competition in the commercial aircraft industry 
has resulted in a renewed interest in high-lift aerodynamics. High-lift configurations for 
commercial transports are characterized by multiple bodies generally in close proximity 
of each other, resulting in complex, highly viscous flows involving merging turbulent 
shear layers, laminar separation bubbles, and extensive regions of separated flow. The 
combination of landing and takeoff flight conditions and aircraft size result in a high 
Reynolds number flow at a relatively low free-stream Mach number. Circulation around 
these configurations can be so great that supersonic flow can exist even at low free- 
stream Mach numbers. 

Given the above characteristics, it is clear that a computational method used to 
simulate such flows must itself have certain features. First, it must be able to deal with 
multiple bodies. This dictates that either unstructured grids or block-structured grids 
be used. Second, it must resolve the flow features, namely the merging shear layers, 
among others. This would require either extreme global grid refinement or an adaptive- 
grid capability. Third, it must take into account compressibility effects. Finally, to 
capture viscous effects, turbulence and transition modelling are required. While either 
grid methodology is capable of performing the above modelings, the focus of this study 
is on an unstructured-grid method. 

Since it is relatively easy to implement an adaptive strategy with unstructured grids, 
these methods have been very popular for solving high-lift flows. Their main drawbacks, 
however, are the memory overhead associated with storing grid connectivity information 
and the computer time associated with indirect addressing. In addition, due to the lack of 
grid structure, it is difficult to implement simple implicit schemes such as approximate 
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factorization, while explicit schemes suffer from slow convergence. Present solvers are 
either implicit schemes utilizing iterative matrix solvers[l][2] or explicit schemes using 
acceleration techniques such as implicit residual smoothing and multigrid[3 ] [4] . 

The multigrid method is particularly appealing since, in theory, the number of 
iterations required to reach a given level of convergence can be independent of the 
grid size. In other words, the work required to achieve a given level of convergence 
depends linearly on the number of grid points. Past efforts with structured-grid solvers 
have shown that remarkable gains in efficiency can be achieved through the use of a 
multigrid algorithm[5][6]. However, the implementation of a multigrid algorithm in an 
unstructured-grid environment is much more difficult. The lack of directionality and 
structure in the grid makes grid coarsening somewhat ambiguous, and intergrid transfers 
are not straightforward. 

The recent effort of Mavriplis[3] to use multigrid with an unstructured-grid solver 
has been very successful. In his method, the solver and grid generator are closely 
coupled. Structured grids are generated around each solid body in the flow field, and the 
resulting points are overlaid and triangulated. For viscous meshes, the structured grid 
is maintained near the surface so that interpolation coefficients can be calculated. The 
scheme is a multistage Runge-Kutta scheme with residual smoothing and the various 
levels of grid refinement used in the multigrid algorithm are generated independently. 

The work presented here is an implementation of multigrid acceleration for an existing 
implicit upwind solver using several of Mavriplis’ techniques. Modifications have been 
made to eliminate the need for grid structure near the surface and hence uncouple the 
solver from the grid generation process. The details of this work are presented in the 
following chapters. 

First is a more detailed description of multigrid methods as well as details of the 
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implementation used in this work. Next, a model problem is described that was used to 
test the multigrid method. Following this is a description of the original implicit upwind 
flow solver, followed by results showing the marked improvement in efficiency obtained 
by using multigrid. 
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n. MULTIGRID 


Many relaxation schemes damp high-frequency error components relatively quickly, 
but are generally slower to damp low-frequency components. By interpolating the 
solution to a coarser grid, these low frequency errors appear as higher frequencies that 
can be damped well by the relaxation scheme. The coarse grid can be used to compute 
a correction to the fine-grid solution to eliminate its low-frequency errors. By using 
successively coarser grids recursively, lower and lower frequency components of the 
fine-grid solution error can be eliminated, and by performing a direct solution on the 
coarsest grid, the convergence rate of the multigrid cycle can be the same as that of the 
relaxation scheme for only the high frequencies. 

In the following sections are descriptions of the basic multigrid methods for both 
linear and nonlinear equations, details of common multigrid cycles and their implemen- 
tation using recursion, and specific details on the implementation of the intergrid transfer 
operators used in the current work. 

Linear Systems 

A system of linear equations can be written as 

L(u) = f (1) 

where L is a linear operator, u is the solution vector, and / is a forcing function. The 
discrete approximation of the system on a grid characterized by spacing h is written 

!*(«*)=/* ( 2 ) 

where u h is the exact solution to the discrete system. Let u h be the current approximation 
to the exact solution u ^ and now define the error v h as 



Now equation 2 can be written as 


L h (u h + v h> j=f h (4) 

which, since L is a linear operator, can be rewritten as 

L h (u h ')+L h (v h ')=f h (5) 

The error v h can be represented on a coarser grid characterized by spacing 2 h provided 
that it is sufficiently smooth to prevent aliasing of high-frequency components on the 
coarse grid. 

An approximation to v h can be calculated on the coarse grid by writing equation 5 
for the coarse grid 

L 2h v 2h = (6) 

where iV 1 is referred to as the restriction operator, which transfers quantities from the 
fine grid h to the coarse grid 2 h. The implementation of this operator is described in 
a following section. Equation 6 can be simplified by defining f 2h = I 2h (f h — L h u h ) 
to obtain 

L 2h v 2h = f 2h (7) 

Once v 2h is obtained, the fine grid solution can be corrected using 

U) = (u h ) + I h 2h v 2h (8) 

where I% h transfers quantities from the coarse grid to the fine grid and is called the 
prolongation operator. Details of this operator are presented in a later section. 

Low-frequency error components can be efficiently eliminated on coarse grids at a 
fraction of the cost of a fine grid calculation. Eliminating these error components on the 
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fine grid is very costly, as many more relaxation cycles are required than would be on 
the coarse grid. In addition, this process can be performed recursively on successively 
coarser grids (i.e. with spacings Ah, 8 h, etc.) with each coarse grid used to compute a 
correction to the next higher grid level. Details on this recursive process are presented 
in the section titled “Multigrid Cycles.” 

Nonlinear Systems 

For systems of nonlinear equations, the step taken between equations 4 and 5 in the 
previous section cannot be performed, so a different formulation must be used. Following 
is a description of the Full Approximation Storage (FAS) scheme [7]. 

Starting with equation 2, subtract L h u h from both sides to obtain 
L h (u h + - L h u h = f h - L h u h = R h 

Written for the coarse grid, this equation becomes 

L 2h ( ll h u h + w 2 ' 1 ) - L 2h ( ll h u h ) = /f (/* - L h u h ) 

By rearranging terms and defining the coarse grid forcing function as 

/“ = /?(/'•- iV) + i“(/2V) 

equation 10 can be written as 

L 2h u 2h _ f2h 

Once u 2h is calculated, the fine-grid solution is updated using 

Note that the difference in formulations does not preclude the use of the recursive 
processes referred to in the preceding section and described in the next. 


(9) 

( 10 ) 

01 ) 

( 12 ) 

(13) 
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Multigrid Cycles 


The recursive formulations referred to in the preceding sections are described below 
and closely follow those of Briggs [8], A particular implementation of a recursive coarse- 
grid correction scheme is referred to as a multigrid cycle. Two specific multigrid cycles 
used in the present work are described along with a generalization of those cycles into 
a single procedure. 

The simplest multigrid cycle is one involving only two levels. A schematic diagram 
of such a cycle is illustrated in figure 1 . Each cycle begins with one or more relaxation 
sweeps on the fine grid. Next, the restriction operator transfers the residual and solution 
vectors to the coarse grid. One or more relaxation sweeps are then performed on the 
coarse grid. Finally, a correction is prolonged from the coarse grid to the fine grid and 
applied to the fine grid solution. This cycle works most efficiently when the coarse grid 
is solved directly, but for most cases, this is still impractical. 

The two-level cycle can approach its maximum performance without a direct solution 
on the coarse grid by using what is known as a V cycle. The relaxation step on the 
coarse grid is now replaced by another two-level cycle. This repeats recursively until 
the coarsest grid available is reached. An example of the resulting cycle for four grid 
levels is illustrated in figure 2 . During the first half of the cycle at each intermediate 
grid level, residual and solution vectors are received from the finer grid, the solution 
is relaxed a given number of times, and residual and solution vectors are then passed 
to the coarser grid below. When the coarsest level is reached, a correction is passed 
successively upward until finally reaching the finest grid. 

A further improvement to the two-level cycle can be made by replacing the coarse- 
grid relaxation with a pair of two-level cycles. This is again done successively at the 
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coarser levels. The resulting cycle is called a W cycle and an example having four grid 
levels is illustrated in figure 3. 

The two cycles above can be combined into a single recursive procedure by specifying 
a cycle index, /t. The resulting generalized cycle is referred to as a fi cycle. Since the 
ultimate result of a multigrid cycle is a correction to the fine grid, the cycle can be 
expressed as a function whose parameters are the current residual and solution vectors 
and whose result is the new solution vector. Stated mathematically, 

u h (14) 

The f.t cycle can now be described by the following recursive procedure: 

1. Relax n times on grid h. 

2. If grid h is not the coarsest level, 

a. Restrict residual and solution vectors to grid 2 h. 

b. Evaluate u 2h <— /.i 2h (R 2h , u 2h ) /i times. 

c. Prolong correction from grid 2h to grid h. 

Appendix A gives details of the coding of this procedure. 

Intergrid Transfers 

The restriction of the solution from a fine grid to a coarser grid and the prolongation 
of the correction from a coarse grid to a finer grid both utilize bilinear interpolation. 
The procedure is described here for the restriction operator, but the process is identical 
for prolongation. 

To get information to the coarse grid, a bilinear interpolation is performed using the 
data at the three vertices of the fine grid cell that encloses each coarse grid node. The 
coordinates of the vertices along with the quantity being interpolated form a plane in 
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a three-dimensional space. Finding the value of the transferred quantity at the coarse 
grid node amounts to solving the equation of the plane at that node. This plane can be 
expressed mathematically as: 

Ax + By + C = q (15) 


where A, B , C are constant coefficients and q is the quantity being transferred. The 
coefficients are determined by assembling a system of equations using the data at the 
cell vertices as follows: 

’*1 2/1 11 Ml [tfi 

X 2 2/2 1 B = <72 ( 16 ) 

.*3 2/3 lj LC'J L?3. 

where the subscript denotes a particular vertex of the fine grid cell. Solving this system 
yields the following expressions for the coefficients: 

<7i(2/2 - 2/3) + <72(2/3 ~ 2/i) + <73(2/1 ~ 2/2) 


A = 


B = 


C = 


21(2/2 - 2/3) + 22(2/3 - 2 /i) + 23(1/1 - 3/2) 

£1(23 - £2) + 92(21 - 23) + <73(3:2 - 21) 

21(2/2 - 2/3) + 22(2/3 - 2/1) + 23(2/1 ~ 2/2) 

91(222/3 - 232/2) + 92(232/1 ~ 212/3) + 93(212/2 ~ 222/1) 


(17) 

(18) 
(19) 


21(2/2 - 2/3) + 22(2/3 - 2 /i) + 23(1/1 - 1/2) 

When equation 15 is evaluated using the coordinates of the coarse grid nodes, it is 
apparent that the value of the quantity at the coarse grid node may be written as the 
sum of geometric weights multiplied by the values of the quantity at the vertices of the 
enclosing fine grid cell, i.e. 


9n = Wiqi + W 2 q 2 + W393 


( 20 ) 


Given the coordinates ( x n ,y n ) of the coarse grid node and the coefficients A,B,C 
substituted into equation 15, the quantities W\,W 2 , W 3 are found by inspection to be 

2„( 1/2 - 2/3) + 2/n(23 “ 22 ) + (222/3 ~ 232/2) 


Wi = 


2l(2/2 - 2/3) + 22(2/3 - 2/l) + 23(2/1 - 2/2) 


( 21 ) 
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( 22 ) 


Xnjm - yi) + Vnjx 1 - £ 3 ) + (£32/1 ~ £l2/3) 
£l('2/2 “ 2/3) + ^ 2 ( 2/3 - 2/l) + £3(2/1 - 2/2) ’ 

£ n (2/1 ~ 2/2) + 2/n(£2 ~ £l) + (£12/2 ~ £22/l) 
£l(2/2 - 2/3) + £2(2/3 ~ 2/l) + £3(2/1 - 2/2) 


(23) 


It is easily verified that the sum of these three weights is unity. 

If the linear interpolation outlined above is used for the transfer of residuals from a 
fine grid to a coarser grid, a situation may arise as shown in figure 4, where a nonzero 
residual at fine grid node P is not utilized on the coarse grid, since none of the fine-grid 
cells having node P as a vertex enclose any coarse-grid nodes; hence, much of the benefit 
of multigrid is lost. In addition, the residual is actually the surface integral of the fluxes 
around the boundary of the control volume and is therefore related to the time rate of 
change of the conserved variables. In order for this rate of change to be the same for 
all grids, it is necessary that the residual transfer be conservative, that is, that the sums 
of the residuals on the fine and coarse grids be equal. For these reasons, the restriction 
process for residuals is handled in the following manner. 

For a given fine grid node, the coarse grid cell that surrounds the node is determined. 
The residual for the fine grid node is then distributed to the vertices of the surrounding 
coarse grid cell. The weights used for the distribution are the same weights used in the 
linear interpolation from the coarse grid to the fine grid. This process ensures that all 
fine-grid residuals contribute to the coarse grid, and that the total residual is conserved, 
since the weights multiplying the residual at any given fine-grid node sum to unity. 

To transfer information from one grid to another using the above interpolation 
operators requires knowledge of which cell of one grid encloses each node of the other 
grid. To determine this information, a tree search similar to that used by Mavriplis[9] is 
used. In this procedure, the nodes of the first grid are first ordered into a list such that 
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a given node has an immediate neighbor appearing earlier in the list. The search then 
proceeds as described in the following paragraph. 

For the first node, an arbitrary cell of the second grid is chosen to start the search. 
If the cell does not enclose the node, the immediate neighbors of the cell are added to 
a list of cells to check (provided the neighboring cells have not already been checked). 
Next, the neighbors of the neighboring cells are checked, and so on until the enclosing 
cell is found. For the remaining nodes, the cell enclosing a neighboring node appearing 
earlier in the list of nodes (i.e. one whose enclosing cell is known) is used as a starting 
point for the search. 

The search algorithm can encounter problems near boundaries, where the grid is 
actually a planar discretization of a curved surface, as illustrated in figure 5. Alternating 
nodes on the fine-grid boundary are displaced away from the coarse-grid boundary. The 
problem is aggravated by concave surfaces and highly stretched viscous meshes, where 
several interior nodes may lie outside the coarse-grid domain, and fine-grid interior nodes 
very close to the surface may receive interpolated information from coarse-grid cells 
farther away from the surface. This latter case is illustrated in figure 6. 

Simple structured-grid algorithms perform intergrid transfers in computational space, 
where grid lines match nicely and operators are straightforward. The equivalent situation 
in physical space is illustrated in figure 7. The following procedure is a way of 
approximating this behavior on unstructured grids by preserving the distance to the 
boundary for each node in a prescribed region. The procedure is described for a fine- to 
coarse-grid transfer. The reverse operation is similar. 

First, the list of boundary faces for each grid is sorted such that adjacent faces are in 
order and the face with the greatest x -coordinate is first in the list. Then, starting at the 
first boundary face of each grid, the boundaries are matched by determining which coarse- 
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grid face is closest to each fine-grid node. Each fine-grid boundary node is assigned 
interpolation coefficients by projecting the node onto the coarse-grid face and computing 
a linear interpolation along the face. The physical displacement required to move the 
fine-grid node to the coarse-grid face is also stored for later use. 

Next, a region near the fine-grid surface is defined in which nodes will be shifted to 
maintain their position relative to the boundary. This is done by first tagging the nodes 
of the grid lying on viscous boundaries. The edges of the grid are then cycled through 
a prescribed number of times. Each cycle through the edges, nodes neighboring tagged 
nodes are themselves tagged. The result for the particular case of a triangulated structured 
grid is that a certain number of layers of grid points have been tagged. In general, nodes 
in a region surrounding the viscous boundary nodes will be tagged. 

Each tagged node is then associated with the boundary face it is nearest. The 
node is projected onto the boundary face, and the previously computed boundary-node 
displacements are then used to determine the displacement to be applied to the interior 
node via linear interpolation. Note that these shifted node coordinates are used only in 
the above procedure for calculating the interpolation coefficients and not in the rest of 
the flow calculation. 
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III. THE MODEL PROBLEM — LAPLACE’S EQUATION 


For the purposes of developing and debugging code for the multigrid cycle and 
its associated intergrid transfer mechanisms, it is beneficial to decouple the difficulties 
associated with the numerics of a problem involving a system of nonlinear equations 
from those associated strictly with information transfer between grids. To this end, 
the two-dimensional Laplace equation (V 2 4> = 0) is coded and used as a test vehicle. 
The boundary conditions are chosen such that the problem represents two-dimensional 
nonlifting potential flow. Namely, ^ = 0 on the airfoil surface(s) and <f> = x at the 
outer boundary (i.e. free stream conditions). 

Spatial Discretization 

The equation is solved in integral form, i.e. 

J V 2 <t>dA = 0 (24) 

s 

At each node, the equation is integrated over a control volume consisting of the cells 
surrounding the node. This control volume surrounding a given node 0 is shown as the 
dark outline in figure 8 and is known as the median-dual control volume. The median 
dual mesh around a given node is formed by joining the centroids of the cells surrounding 
the node and the midpoints of the edges joining to the node. Using Green’s theorem, 
the integral equation becomes 

£v<f>-hdl = 0 (25) 

C 

and the integral over the control volume can be evaluated by integrating around its 
boundary. 

In reference 1, it is shown that integrating over the control volume described above is 
equivalent to a Galerkin finite-element approximation. Now contributions from each edge 
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that connects to a given node can be considered individually given the nodes that form 
the edge and the centroids of the cells on either side of the edge. Each edge contributes 
two segments to the boundary of the median-dual control volume, as shown in figure 9. 

The Laplacian operator can now be discretized as the sum of weights associated with 
edges of the mesh multiplied by the difference in the solution along the edge, i.e. 

( V ^)o = £ M ( 26 ) 

i&Xo 

where xo is the set of edges connecting to node 0 and Mi is the edge weight. Reference 
1 presents a detailed derivation of these edge weights for the Laplacian operator, and a 
brief summary of the results is presented here. Using the notation of figure 9, the edge 
weights can be expressed as: 

1 ^j-i-i/2 ' ii'i+i 1/2 • *W-i 

Mi = -- 7 

4 [ A i + 1/2 1/2 

where ^ 4 ,_i /2 and ^4 ,+i/ 2 are the areas of the cells to the right and left of the edge 
connecting node 0 and node i, respectively. Identities are then used to express this in 
terms of the single edge. The resulting expression is: 

Note that these weights depend only on the geometry and not on the solution, so they 
can be precomputed and stored for the duration of the calculation. 

Iteration Scheme 

Jacobi iteration with a relaxation parameter is used to advance the solution. Starting 
with equation 26 , the residual at iteration n is defined as 

RS = Ys - $>) ( 29 > 


1 {fij • hj - 9 • nfl, \hj • n,- - 9n Lt • h Li 

4 A i+ 1/2 "4.-1/2 
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Laplace’s equation is now discretized using <f > o at iteration n + 1 and 4>i at iteration n 
as follows: 


£A/ j (^-C H )=0 

i€\ o 


This equation is then solved for <t> q + 1 to yield 


E 


±n+l _ ie *o 

<Po ” 


E Mi 

i€\ o 

Now subtract </>q from both sides to yield: 

E 

jLTI+1 X 71 

h EW 

i€*o 

The increment to the solution A<^q is then calculated by: 


? = A? + 1 - t hH = 


&<t>0 = K 


as 

E Mi 

*6 Yo 


(30) 


(31) 


(32) 


(33) 


Over- or underrelaxation is accomplished simply by adding a relaxation parameter 
a as follows: 

(34) 

Yo 

If a is greater than unity, the solution is overrelaxed, and if it is less than unity (but 
greater than zero), the solution is underrelaxed. Of course, Jacobi iteration is recovered 
when ce is equal to unity. 

Alternately, a red-black scheme may be used. In this scheme, the grid is divided into 
two “colors” — red and black — depending on whether the node number is even or odd, 
respectively. The increment to the solution is calculated as follows: 

1. Calculate the residual for all nodes. 
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2. Update nodes colored red. 

3. Recalculate the residual for all nodes. 

4. Update nodes colored black. 

The scheme requires an additional residual calculation at each iteration, but exhibits better 
smoothing properties from the perspective of multigrid methods. [8] 

Results 

The Laplace solver was run in different modes for two cases. The first case is the 
simple case of a square with homogeneous Dirichlet boundary conditions on all its edges. 
The grids consisted of Cartesian meshes with each cell cut across one of its diagonals 
to form triangles. An example of one of these meshes is shown in figure 10. An initial 
condition of unit potential was used on the interior grid. Finally, the solver was used to 
simulate nonlifting potential flow over a NACA 0012 airfoil at zero angle of attack. 

Square Domain 

The first study, shown in figure 1 1, is a demonstration of the effect of the relaxation 
scheme on the performance of the multigrid scheme. All three cases were run with a 
direct solution on the coarsest grid and used a 4-level V cycle. Clearly, the Jacobi scheme 
is inferior, as it poorly damps the high frequency components of the solution error. The 
damped Jacobi scheme (a=0.5) exhibits better performance since underrelaxation greatly 
improves the damping of the high frequencies[8]. The red-black scheme shows excellent 
performance. The damped Jacobi scheme is used for the remaining studies on the square 
domain. 

In practice, a direct solution on the coarsest grid is not used, but by performing more 
relaxations on the coarse grid, improving the cycle, and using additional grid levels, the 
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performance of the scheme can approach that which is achieved by a direct solution on 
the coarse grid. The following studies illustrate this point. 

The next study shows the effect of the level of convergence of the coarsest grid on 
the convergence rate of the multigrid scheme. Figure 12 shows convergence histories for 
a 4-level V cycle. The parameter MITER is the number of damped Jacobi relaxation 
sweeps performed on the coarsest mesh during each multigrid cycle. The curve marked 
MITER -*• oo corresponds to a direct solution on the coarsest grid. It is obvious that a 
relatively meager increase in the number of relaxation sweeps on the coarsest grid has a 
profound impact on the convergence of the multigrid scheme. Although not shown, the 
convergence history for M ITER = 10 is virtually identical to that of a direct solution 
on the coarsest mesh. These extra sweeps on the coarsest mesh are cheap in terms of 
computational work in comparison with the finest mesh. It will be shown, however, that 
the same increase in performance can be achieved by changing the multigrid cycle. 

Figure 13 shows that using the W cycle instead of the V cycle results in the same 
improved convergence of the multigrid cycle. Both cases use 4 grid levels with 1 
damped-Jacobi relaxation sweep on the coarsest mesh. The residual is now plotted 
against computer time to reflect the actual computational work done, since a W cycle 
requires more computational work per multigrid cycle than does a V cycle. 

In figure 14, the benefit of using additional grid levels in the multigrid cycle is 
apparent. All cases use a V cycle with 10 damped-Jacobi relaxation sweeps on the 
coarsest grid. Each level of coarsening allows the relaxation scheme to damp lower and 
lower frequency components of the solution error more effectively. This figure along 
with figure 12 shows that the scheme can approach its best performance without a direct 
solution on the coarsest grid. 
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NACA 0012 Airfoil 


To demonstrate the ability of the Laplace solver to calculate nonlifting incompressible 
flows, results for a NACA 0012 airfoil at zero incidence are presented. The red-black 
scheme has been used for this case with one relaxation sweep at each grid level. At the 
airfoil surface, the derivative of the potential normal to the surface vanishes, and free 
stream potential is enforced at the far field boundary. Four levels of grid refinement were 
generated, the finest consisting of 14,269 nodes. A portion of this grid is shown in figure 
15. Table 1 summarizes the grid levels used for these cases. 


Table 1. — Summary of grid sizes for grids around NACA 0012 airfoil for inviscid flow. 


Grid level 

Total nodes 

Nodes on surface 

0 

14,269 

256 

l 

3796 

128 

2 

1081 

65 

3 

424 

36 


Figure 16 shows convergence histories with and without multigrid. The multigrid 
algorithm has clearly provided a substantial improvement in the rate of convergence. 

A comparison of surface pressure coefficient is presented in figure 17 along with an 
analytical solution[10]. The solution agrees well with the analytical data. 
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IV. VISCOUS FLOW SOLVER 


The solver used in the present work was developed by W. Kyle Anderson at the 
NASA Langley Research Center. A summary of the salient features of the code are 
presented in the following chapter for completeness. A more detailed description of the 
solver is presented in appendix B, and further information can be found in reference 1 1. 

The governing equations are the time-dependent two-dimensional Reynolds-averaged 
Navier-Stokes equations in conservation-law form, which are integrated in time to obtain 
a steady-state solution. Ideal-gas assumptions are made, and either of the one-equation 
turbulence models of Spalart and Allmaras[12] or Baldwin and Barth[13] may be used 
for calculating turbulent flows. 

The temporal formulation consists of a backward-Euler time difference, with the 
resulting linear system of discrete equations being solved iteratively using red-black 
relaxation. The result is that at each iteration of the nonlinear system, a prescribed 
number of “subiterations” are performed to obtain an approximate solution to the linear 
system. A finite-volume formulation is used to discretize the governing equations in space 
at each node. For the convective and pressure terms, the upwind scheme of Roe[14] is 
used, while a simple central difference is used for the viscous terms. 

The same scheme is used to solve the turbulence equation; however, this calculation 
is carried out separately. At each iteration, a prescribed number of subiterations are 
performed on the flow equations while holding the turbulence fixed, followed by an 
update of the turbulence equation holding the flow quantities fixed. 
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V. RESULTS 


The following chapter presents results in the form of histories of the temporal 
convergence of both the L 2 norm of the residual for the continuity equation and the 
lift coefficient, C/. These quantities are shown versus computer time. Unless otherwise 
noted, all cases were run on a Cray-YMP, and grids were generated using the advancing- 
front method described in reference 15. 

First, cases of inviscid and laminar flow over a NACA 0012 airfoil are presented, 
followed by several cases of turbulent flow over both the RAE 2822 airfoil and a 
3-element airfoil. Finally, results are presented for a turbulent case on a grid for which a 
calculation without multigrid is impractical. All turbulent cases used the Spalart-Allmaras 
turbulence model. 

Euler Solution 

The first case presented is that of inviscid, transonic flow over a NACA 0012 airfoil. 
The free stream Mach number is 0.8, the angle of attack is 1.25°, and the same grids 
that were used for the potential-flow case presented earlier are used (see table 1). Figure 
18 shows the convergence histories for several V cycles in comparison to the original 
scheme. Note that the multigrid scheme substantially improves the convergence rate, and 
that the improvement increases as more grid levels are used. Results for the W cycle 
are shown in figure 19. Note again the substantial improvement in convergence. Figure 
20 shows the best V-cycle result and the best W-cycle result together with the results 
for the scheme without multigrid. The V cycle and W cycle perform similarly versus 
computer time; however, the relaxation parameters used were those found to work well 
for the base scheme (i.e. without multigrid), specifically, 20 subiterations were used. 
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If the scheme is used without multigrid, the tunable parameters (i.e. CFL number and 
number of subiterations) must be chosen to give the fastest convergence. With a multigrid 
method, only the high-frequency error components need to be damped quickly on all but 
the coarsest grid. This would seem to imply that with the present scheme, a further 
reduction in computer time could be achieved by reducing the number of subiterations. 
Figure 21 shows the effect of reducing the number of subiterations for a four-level 
W cycle. While convergence per cycle is slightly compromised, convergence versus 
computer time is improved due to the decrease in computational work per cycle. 

Residual and lift histories for the V and W cycles using five subiterations at each 
grid level are shown in figure 22. In this case, the W cycle slightly outperforms the V 
cycle, while both obtain a steady value for the lift approximately four times faster than 
the base scheme alone. 

Laminar Navier-Stokes Solution 

Figure 23 shows convergence histories for a case of laminar flow over a NACA 0012 
airfoil at an angle of attack of 3°, a free-stream Mach number of 0.5, and a Reynolds 
number of 5000. The grids used in the multigrid cycle are shown in figure 24 and are 
summarized in table 2. The two multigrid cycles achieve a steady lift coefficient in a 
Table 2. — Summary of grid sizes for grids around a NACA 0012 airfoil for laminar flow. 


Grid level 

Total nodes 

Nodes on surface 

0 

16,116 

256 

1 

5004 

128 

2 

1891 

65 

3 

1237 

40 


fraction of the time taken by the original solver, and the W cycle has a slight edge over 
the V cycle, particularly in the convergence of the lift coefficient. 
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Turbulent Navier-Stokes Solutions 


After presenting a case of transonic flow over an RAE 2822 airfoil, several cases 
of flow past a 3-element airfoil are shown, including a case previously impractical to 
calculate. For all 3-element cases presented, the free-stream Mach number is 0.2 and 
the Reynolds number is 9 million. 

Figure 25 shows convergence histories for flow past an RAE 2822 airfoil at 2.81° 
angle of attack, free-stream Mach number of 0.75, and a Reynolds number of 6.2 million. 
The grids were generated using the method described in reference 16, and a summary 
of their characteristics is presented in table 3. The residual for both multigrid cases 
Table 3. — Summary of grid sizes for grids around RAE 2822 airfoil for turbulent flow. 


Grid level 

Total nodes 

Nodes on surface 

0 

13,385 

208 

1 

3359 

104 

2 

847 

52 

3 

219 

26 


converges a few orders of magnitude before cycling about a nearly constant level. Other 
runs have shown that this phenomena is a result of an adverse coupling of multigrid 
and the turbulence model, as holding the turbulence quantity constant after some level 
of convergence has been reached causes the residual to continue decreasing. Steady 
lift for both multigrid cases is still achieved prior to the cycling of the residual and in 
significantly less time than for the original scheme. 

Figures 26 and 27 show the distributions of the surface pressure coefficient and 
skin friction coefficient, respectively, for the 4-level W cycle along with experimental 
data[17]. The computed results are in good agreement with the experimental data, and 
are virtually identical to results obtained with the base scheme. 
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A case of turbulent flow over the 3-element airfoil shown in figure 28 at 16.21° 
is shown in figure 29, and the characteristics of the grids are summarized in table 4. 
The multigrid cases again cycle about some level after a certain level of convergence is 


Table 4. — Summary of grid sizes for grids around a 3-element airfoil for turbulent flow. 


Grid level 

Total nodes 

Nodes on surfaces 

0 

97,088 

1340 

1 

34,987 

671 

2 

14,278 

340 

3 

6657 

178 


reached. The lift, however, converges for all three cases and does so much more rapidly 
for the multigrid cases, with the W cycle having a significant edge. 

The same configuration at a higher angle of attack is shown in figure 30. The precise 
angle of attack is 21.34° and is near maximum lift as determined by experiment! 18], 
The spikes in the residual histories are a result of restarting the code. Specifically, a 
point-vortex is applied at the outer boundary whose strength depends on the lift, which 
is not available during the first iteration since it is presently calculated after the residual. 
This is easily cured by computing the lift before computing the residual. Note that for 
this run, the V-cycle case continues converging while the W-cycle residual again cycles 
after less than two orders of magnitude of convergence. The multigrid scheme again 
shows considerable improvement over the base scheme. 

To further demonstrate the advantages of multigrid, the 3-element airfoil was run at 
an angle of attack of 16.21° on a grid consisting of 309,770 nodes. The characteristics 
of the full set of grids is given in table 5. This case had been considered impractical 
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Table 5. — Summary of grid sizes for grids around a 3-element airfoil for turbulent flow. 


Grid level 

Total nodes 

Nodes on surfaces 

0 

309,770 

2679 

1 

97,088 

1340 

2 

34,987 

671 

3 

14,278 

340 


with the original solver due to the nonlinear increase in computer time required with the 
increase in the number of grid points. The convergence histories are shown in figure 
31. Computer restrictions dictated that only 50 cycles could be calculated in a single 
run. As explained earlier, the spikes in the convergence histories are a result of restarting 
the code. The W cycle exhibits oscillatory behavior in the residual, while the V cycle 
continues converging. The lift seems nearly steady for the W cycle, but when viewed on 
a smaller scale, it exhibits small-scale oscillations. The lift for the V cycle, however, is 
steady, and the surface pressure distributions for this case are presented in figure 32. 
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VI. CONCLUSIONS 


A multigrid algorithm has been implemented in an existing code for solving turbulent 
flows on triangular meshes. Intergrid transfer operators have been used that ensure 
conservation of the residual and preserve smoothness of the solution near solid surfaces. 
Once coded, the multigrid algorithm and intergrid transfer operators were used to solve 
Laplace’s equation to verify correct operation. 

The Laplace solver with the red-black relaxation scheme and multigrid algorithm 
is very efficient for solving nonlifting potential flow on unstructured grids, and was 
indispensable for validating intergrid transfer operators and the multigrid cycle itself. 

The multigrid algorithm has improved convergence significantly for both inviscid and 
laminar viscous flows. For the turbulent flows, the improvement with multigrid can be 
quite dramatic, with increasing improvement with grid refinement. 

Several avenues of future study exist as a result of this work. The apparent adverse 
interaction between the W cycle and the turbulence model will require a significant effort 
to resolve. The method can also be extended to three-dimensions, or to higher-order 
methods. 
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APPENDIX A: CODING OF MULTIGRID CYCLE USING RECURSION 


Since the main driver of the flow solver is written in C, which allows a function 

to call itself recursively, implementation of the /t cycle described in the text is very 

straightforward, and can be nearly literally translated into C-code. Following is the code 

fragment representing the multigrid cycle: 

mu eye (mu , i f ine , igridl , igrid2 , grid, miter ) 

GRID *grid; 

int mu, if ine, igridl , igrid2 ; 
int *miter; 

{ 

int i , j ; 

relax (miter [ igridl] , grid[igridl] ) ; 
if (igridl == ifine) 

{ 

f 77L2NORM ( ) ; 
f 77FORCE ( ) ; 

} 

if (igridl < igrid2) 

{ 

restricter (grid [igridl] , grid[igridl+l] ) ; 

for (i = 0; i < mu; ++i) 

{ 

mucyc (mu, ifine, igridl+1 , igrid2 , grid, miter) ; 

} 

prolong (grid [igridl] , grid[igridl+l] ) ; 

} 


> 

In this routine, f 77L2NORM and f77FORCE are FORTRAN routines that calculate 
quantities used to monitor convergence on the fine grid. These two routines have many 
arguments, but they are omitted here for clarity. Note that grids are denoted by index 
numbers (0, 1, 2, ...) rather than characteristic spacings (h, 2h, 4 h, ...). The parameter mu 
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is the cycle index, while igridl and igrid2 are the finest and coarsest grid levels in 
the cycle, respectively. The parameter if ine is a copy of the initial value of igridl. 
The argument grid is an array of structures having one entry for each grid level. Each 
structure contains parameters indicating the size of the corresponding grid, as well as 
pointers to arrays containing connectivity and field information. The argument miter is 
an array containing the number of relaxation sweeps to be performed at each grid level. 
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APPENDIX B: DESCRIPTION OF ORIGINAL VISCOUS FLOW SOLVER 


Governing Equations 

The relaxation scheme solves the Reynolds-Averaged Navier-Stokes (RANS) equa- 
tions in conservation-law form. These equations are given in vector form by 

A m + f Fj • hdl - j F v • ndl = 0 (35) 

an an 

where n is the outward-pointing unit normal to the surface of the control volume 3D. Q 

is the vector of conserved state variables given by 

P 

I 

Q = 


pu 

pv 

E 


(36) 


and F, and F v are the inviscid and viscous fluxes, respectively, through the surface of 
the control volume dfl and are given by 


Fi = fi+gj = 
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The shear stress and heat conduction terms in the viscous fluxes are given by 
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The perfect-gas equation of state is used to define the pressure p and is given by 


P =(7-1)[£-p(« 2 +» 2 )/2] 


( 44 ) 


and the laminar viscosity p is given by Sutherland’s law 

P (1 "b C* ) /'ti /rp \3/2 


Poo ( T/T, x + C *) 


{T/TooY 


(45) 


where C* = is Sutherland’s constant divided by a free-stream reference temperature 
assumed to be 460° Rankine. 

The eddy viscosity p^ is obtained by either of two one-equation turbulence closure 
models. The first, developed by Baldwin and Barth[13], is derived from the k-e equations. 
The second, developed by Spalart and Allmaras[12], relies more heavily on empiricism 
and dimensional analysis. The turbulence model is solved separately from the rest of the 
system, but uses the same solution scheme, and, although multigrid is also used with the 
turbulence model, it remains decoupled from the rest of the system. The Spalart-Allmaras 
model is used for all turbulent calculations in this study. 


Time Integration 

The governing equations are integrated in time to the steady-state solution using 
a linearized backward-Euler time-differencing scheme. The resulting system of linear 
equations can be expressed as 

[A] n {AQ} n = {R}” ( 46 ) 


where 


4 3R” 

[A1 = Tt 1 + -gQ 


(47) 


The solution of this linear system is obtained iteratively via a classic relaxation procedure. 
To differentiate between the nonlinear and linear systems, the term “iteration” is used to 
refer to the nonlinear system, while “subiteration” is used to refer to the linear system. 
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To illustrate the scheme used, let the matrix [A] n be written as the sum of two 
matrices representing the diagonal and off-diagonal terms 

[A] n = [D] 71 + [Of (48) 

The simplest method for solving the linear system is commonly referred to as Jacobi iter- 
ation and consists of moving all off-diagonal terms to the right-hand side and evaluating 
them at the previous subiteration i. This can be written as 

[D] n {AQ} 1+1 = [{R} n - [0] n {AQ}'j (49) 


The convergence of this method is accelerated somewhat using a red-black scheme 
where even-numbered nodes are update using the Jacobi scheme just described, followed 
by the odd numbered nodes using the update values at the even-numbered nodes. This 
scheme can be written as 


[Df{AQ} i+1 = \{R} n - [0] n {AQ}‘*’' 


(50) 


where {AQ} ^ is the most recent value of {AQ} and will be at subiteration i + 1 for 
the even-numbered nodes and at subiteration i for the odd-numbered nodes. 


To further accelerate convergence, local time-stepping is used. A separate time step 
is calculated at each node using the inviscid stability limit. 


Spatial Formulation 

The spatial discretization is a finite-volume formulation in which the inviscid and 
viscous fluxes are integrated over the median-dual control volume surrounding each node 
(see figure 8). Green’s theorem is used to change the volume integrals to surface integrals 
over the edges of the dual mesh. These surface integrals can be calculated using edge 
formulas as described in reference 1. 
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The inviscid fluxes, F;, are obtained on the edges of the control volume using Roe’s 
approximate Riemann solver[14]. The viscous fluxes, F v , are computed using a simple 
central difference. 
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FIGURES 



I - Relaxation sweep(s) 
R - Restriction 
P - Prolongation 


Figure 1: Schematic of two-level multigrid cycle. 



I - Relaxation sweep(s) 
R - Restriction 
P - Prolongation 


cycle for four grid levels. 



Relaxation sweep(s) 

Restriction 

Prolongation 


Figure 3: Schematic of W cycle for four grid levels. 
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Figure 4: Example of a fine-grid node P that will not contribute 
information to the coarse grid if linear interpolation is used. 




Solid surface 
Fine-grid boundary node 
Fine-grid boundary face 
Coarse-grid boundary node 
Coarse-grid boundary face 


Figure 5: Discretization of a curved boundary surface for both a fine and a coarse grid. 


Fine Grid 
Coarse Grid 

Figure 6: Effective interpolation near viscous surfaces. The diagonal edges 
cutting across the quadrilateral cells are omitted for visual clarity. 





Fine Grid 
Coarse Grid 


Figure 7: Effective interpolation near viscous surfaces for structured grids. 
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Figure 8: Median-dual control volume for node 0. 



Figure 9: Contribution of an individual edge 
to the median-dual control volume for node 0. 
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Figure 10: Sample unstructured grid for a square. 



0 20 40 60 80 100 


iteration 


Figure 11: Effect of relaxation scheme on multigrid performance. All runs are 
made using a 4— level V cycle with a direct solution on the coarsest mesh. 




Figure 12: Effect of coarse grid convergence level on convergence of 
the multigrid cycle. Although not shown, convergence for 
MITER = 10 are virtually identical to those for MITER —*■ oo. 



Figure 13: Effect of multigrid cycle on convergence. Both cases use 4 grid 
levels with 1 damped-Jacobi relaxation sweep on the coarsest grid. 
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Figure 14: Effect of number of grid levels on convergence of the multigrid scheme. All 
cases use a V cycle with 10 damped-Jacobi relaxation sweeps on the coarsest grid. 



Figure 15: Portion of grid around a NACA 0012 airfoil. 
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CPU time, sec. 


Figure 16: Convergence history for nonlifting potential 
flow over a NACA 0012 airfoil at zero incidence. 



Figure 17: Surface pressure coefficient distribution 
on a NACA 0012 airfoil at zero incidence. 







Cycle 


CPU time, sec. 


Figure 18: Convergence histories for several V cycles versus both cycle number and 
computer time for inviscid flow over a NACA 0012 airfoil (Moo = 0.8, a = 1.25°). 




CPU time, sec. 


Figure 19: Convergence histories for several V cycles versus both cycle number and 
computer time for inviscid flow over a NACA 0012 airfoil (Moo = 0.8, a = 1.25°). 
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0 100 200 300 



Cycle 


CPU time, sec. 


Figure 20: Comparison of performance of V and W cycles versus both cycle number 
and computer time for inviscid flow over a NACA 0012 airfoil (Mx> = 0.8, a = 1.25°). 




Figure 21: Effect of number of subiterations on performance of a 4-level W 
cycle for inviscid flow over a NACA 0012 airfoil (Moo = 0.8, a = 1.25°). 
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Residual 



CPU time 


CPU time 


Figure 22: Comparison of performance of V and W cycles for 
inviscid flow over a NACA 0012 airfoil (Moo = 0.8, a = 1.25°). 
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Figure 23: Comparison of performance of V and W cycles for laminar 
flow over a NACA 0012 airfoil (Moo = 0.5, a = 3°, Re = 5000). 
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a. Grid level 0. 



b. Grid level 1. 



c. Grid level 2. 

Figure 24: Grids around NACA 0012 airfoil 
used for laminar-flow case. (Continued . . . ) 
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d. Grid level 3. 


Figure 24: Grids around NACA 0012 airfoil used for laminar-flow case. 
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Figure 25: Comparison of performance of V and W cycles for turbulent 
flow over an RAE 2822 airfoil (Moo = 0.75, a — 2.81°, Re = 6.2 x 10 6 ). 


45 





Figure 26: Comparison of surface pressure coefficient distribution with 
experiment on RAE 2822 airfoil (M 0 0 = 0.75, a = 2.81°, Re = 6.2 x 10 6 ). 


xIO " 3 



x/c 

Figure 27: Comparison of skin friction coefficient distribution with 
experiment on RAE 2822 airfoil (M 0 0 = 0.75, a = 2.81°, Re = 6.2 x 10 6 ). 
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Figure 28: Geometry of 3-element airfoil. 
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Figure 29: Comparison of performance of V and W cycles for turbulent flow 

over a 3-element airfoil (M^ = 0.2, a = 16.21°, Re = 9 x 10 6 ). 
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Figure 30: Comparison of performance of V and W cycles for turbulent flow 
over a 3-element airfoil (Moo = 0.2, o; = 21.34°, Re — 9 x 10 6 ). 
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Figure 31: Comparison of performance of V and W cycles for turbulent flow over a 

3-element airfoil (Moo = 0-2, a = 16.21°, Re — 9 x 10 6 ) on a grid of 309,000 nodes. 
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Figure 32: Comparison of distributions of surface pressure coefficients for V 
cycle with experiment for turbulent flow over a 3-element airfoil 
(Moo = 0.2, a = 16.21°, Re = 9 x 10 6 ) on a grid of 309,000 nodes. 
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